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A METHOD OF DETERMINING THE EFFECT OF AIRPLANE 
STABILITY ON THE GUST LOAD FACTOR 
By Bernard Mazelsky and Franklin W. DIederich 

SUMMARY 


The equations of motion of a conventional airplane penetrating a 
gust are given in detail for determining the effects due to stability 
on the gust load factor. A convenient numerical method is derived in 
matrix notation which affords a systematic procedure of solving the 
equations for a unit jump or arbitrary forcing function. The solution 
for the airplane motions for the stick-free condition is modified by 
neglecting the elevator terms so that the effects under similar stability 
conditions may be calculated with the stick in a fixed position. If the 
complete response of the airplane is known, wing and tail loads may be 
computed fairly conveniently. Sample calculations were performed to 
illustrate the application of the equations. 


INTRODUCTION 


In the past the attempts to study the most effective stability 
parameters in a gust have been hampered by the complexity of the equa- 
tions and the time consumed in their solution. In numerous analyses, 
the pitch effects are either considered empirically or neglected 
entirely. A method of calculation that would be feasible, would give 
reasonably accurate results, and would require only a moderate expend! 
ture of time is, therefore, needed for analytical studies of airplane 
motions in gusts. 

Analytical solutions by means of operators have been made for 
calculating the motions of a wing penetrating a sharp-edge gust (see 
reference 1) . An extension of this method for solving the more com- 
plicated equations would be impractical. In reference 2, calculations 
were made by an iterative method for obtaining the solutions of the 
equations of motion in two degrees of freedom of a canard airplane. 

For this analysis the downwash effects could be neglected and the 
equations thereby considerably simplified. 
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The present analysis gives an application of the calculus of 
finite differences for solving the equations of motion of airplanes in 
gusty air. The numerical solutions have been derived in matrix notation 
from the equations of motion in three degrees of freedom - vertical, 
pitch, and elevator displacement - forward-speed variations being 
neglected, and sample calculations are made. For a complete under- 
standing of the method, a working knowledge of matrix operations is 
required. 


SYMBOLS 


L 

CL 

a 

dCL 

da 

A 

P 

U 

V 

q. 

g 

s 

c 

W 


S 1 


lift, pounds 
lift coefficient 
angle of attack, radians 

slope of lift curve, per radian 

aspect ratio 

mass density of air, slugs per cubic foot 
gust velocity, feet per second 
forward velocity, feet per second 

dynamic pressure f 

acceleration due to gravity, feet per second per second 
area of surface, square feet 
reference chord, feet 
airplane weight, pouiids 

mass parameter 

nondimensional distance penetrated into gust, chords 
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s 

As 

in 

r 



$ 


CL 


€ g 


$ 


e 


CLq, 


variable of integration 

increment of s used in matrix solution 

variable in recurrence formula ( 

\As 

nondimens ional discontinuity taken 

+ 1 

lag of tail penetration I — c — s 

normalized unsteady-lift function < 
sharp-edge gust 

normalized unsteady-lift function < 
arbitrary- shape gust ■ 

equivalent normalized unsteady-lift function on tail due to 
downwash caused by wing penetration into a sharp-edge gust 

equivalent normalized unsteady-lift function on tail due to 
downwash caused by wing penetration of an arbitrary- 
shape gust 

normalized unsteady-lift function for a unit jump of angle 
of attack 


to nearest integer due to 


iue to penetration of a 


iue to penetration of an 



e S 

I 


B 

C 

An 

0 


equivalent normalized unsteady-lift function on tail due to 
downwash caused by a unit jump of angle of attack on wing 

normalized unsteady-downwash angle at horizontal tail due to 
penetration of wing into a sharp -edge gust 

normalized unsteady-downwash angle at horizontal tail due 
to penetration of wing into an arbitrary- shape gust 

normalized unsteady-downwash angle at horizontal tail due 
to unit jump of angle of attack on wing 

transformed unsteady-lift function Ci^ for wing ^1 - 

transformed unsteady-lift function for tail ^1 - 

load-factor increment encountered by airplane, multiples of 
acceleration due to gravity 

pitch angle, radians 
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elevator-deflection angle, radians 

horizontal distance from center of- gravity of airplane to 
aerodynamic center of wing, positive when aerodynamic 
center is ahead of center of gravity, feet 

horizontal distance from center of gravity of airplane to 
aerodynamic center of tail, feet 

distance between trailing edge of wing mean aerodynamic 
chord to leading edge of tail mean aerodynamic chord, 
feet 


lag for downwash angle to be effective on tail when wing 
undergoes a unit Jump of~angle of attack 

+ 0.23 


As 

lag for unsteady-lift function due to downwash to be effec- 
tive on tail when wing undergoes a unit jump of angle of 

R* + 0.80> 

attack \ 

\ As 


distance between airplane center of -gravity and elevator 
hinge, feet 


tail-efficiency factor 



steady-downwash angle per unit jump of angle of attack on 
wing 


asymptotic value of downwash angle per-unit jump of angle 
of attack on wing including effects due to unsteady lift 
on tail 

pitching moment of inertia of airplane about center of 
gravity, slug-feet 2 

moment of inertia of elevator about its hinge, slug-feet 2 

moment of inertia of control stick about its pivot, 
slug -feet 2 
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E total hinge moment, positive when moment tends to depress 

trailing edge 


He 

H s 


r 


K 

Q 

C *a 

Ch S 


gravitational moment of elevator about its hinge, positive 
when moment tends to depress trailing edge , slug-feet 

gravitational moment of control stick about its pivot, 
positive when moment tends to depress trailing edge, 
slug-feet 

control gearing ratio, angular stick deflection divided by 
angular elevator deflection 

viscous damping constant, in control system, pounds per foot 
per second 

ratio of damping velocity to elevator velocity, feet per 
second per radian per second 

elevator hinge-moment coefficient due to angle -of -attack 
change on tailj floating tendency is positive when 
surface floats against relative wind 

elevator hinge-moment coefficient due to elevator deflection 
restoring tendency is positive when surface is over- 
balanced 


dn, e 

dS 

(Bi) 

\ZdSI a 




elevator-effectiveness factor 


part of additional lift due to angular velocity of elevator 
caused by acceleration of potential flow 


part of additional lift due to angular velocity of elevator 
caused by effective increase in camber 


part of. hinge moment due to angular velocity of elevator 
caused by acceleration of potential flow 

part of hinge moment due to angular velocity of elevator 
caused by effective increase in camber 


Operational symbols: 


D, 


D 


differential operators 
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single integration of arbitrary function y(s) 

Iiy double integration of arbitrary function y(s) 

Matrix symbols : 

constant elements in matrix solution 

column matrix representing variables in recurrence formula 

column matrix representing the forcing function due to gust 
in recurrence formula 

square matrix representing inverse of the matrix coeffi- 
cients of “ An, lP-6, and D^5 at point m 

. . Js6~] rectangular matrices representing constant elements 
in recurrence formula 

row matrix 

matrices representing constant elements in recurrence 
formula for stick-fixed condition 

Subscripts: 
e elevator 

max maximum 

w wing 

t tail 



DERIVATION OF METHOD 
Differential Equations of Motion 


The evaluation of the effects of stability on the gust load factor 
may be determined from an analysis of the motions of the airplane under 
various stability conditions. In setting up the differential equations, 
the following assumptions are made: 
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(a) The airplane is free to move vertically and to pitch about 
its lateral axis and has an elevator restrained only by the force 
produced by the viscous friction in the control system. 

(b) The airplane is in steady' level flight before entering the 
gust and has initially no displacements, velocities, and accelerations 
with regard to its degrees of freedom. 

(c) The forward speed is constant. 

(d) The forces producing the hinge moment on the elevator are 
steady-state forces. 

(e) The over-all aerodynamic force on the whole tail is a tran- 
sient force. 

An equation of equilibrium exists between inertia and aerodynamic 
forces and moments for each degree of freedom. In figure 1, the 
forces and moments acting on a conventional airplane are shown. All 
distances, forces, and moments are shown in a positive direction. The 
three equations representing the incremental forces and moments become 
(see fig. 1) 

2 

WAn = wi-^i==Lw + Lt=WAn w + W£nt (l) 

dt 2 S 


I ^ jfie = An w - WI t Ant (2) 


Ie (D 2 0 + D%) - He^g An - l h ¥£ p2 0 j 


r fc B ^2 (r D^5 - D 2 S) - E s g An 


= H 


(3) 


The transient aerodynamic forces produced by the airplane 
penetrating the gust and the forces due to its motion relative to an 
equilibrium position may be appropriately described by a set of 
integral equations. That part of. the force due to the penetration of 
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the gust is independent of the motions and is therefore analogous to 
the forcing function acting on any dynamic system. The forces produced 
for each degree of freedom must be accounted for separately. These 
individual forces are described by a number of Duhamel's integrals. 

If superposition is valid, the individual forces may be summed to give 
the total force acting on the body due to its motion. These aero- 
dynamic forces may be converted to their corresponding load-factor 
increments. Accordingly, for the motion of the wing, the load-factor 
increment An w is described by the following integral equation: 


rn tj r s i 

An w (si) = Kp J QLg w ( S 1 - s)D ds - Ki j. J c La w ^ sl " ‘O^ 8 ) ds + 


K 7 f 1 PLo^Csi - s)D 9 ds - Kio J 1 CL^Csp - s)D 2 0 ds (4) 


where the coefficients of the integrals are 


Kl 



| VS 


^max 

W 


(4a) 


*4 = 



(4b) 


k t 


f ^k) SP V 2 
,da / w W 2 


(4c) 


Kio = 



(4d) 


The expression for the forces acting on the tail correspond to 
the terms describing the wing together with the' forces caused by the 
interaction of the wing and the tail. This last effect, called down- 
wash, introduces additional Duhamel's integrals which are functions of 
the motions of the wing. The summations of the integrals represent 
the total aerodynamic forces on the tail and may be converted to load- 
factor increments as follows: 
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P S 1 u r S1 

**(<*■) -KJ o w as + K3 Jo ‘H (sl - 6,D 


u 

^max 


ds - 


X s x _ 

Cto+.Csi- s)An(s) ds - Kg j ~ Ql^si- s)An(s) ds + 


JQ t 


p s i 


si 


.si 


% j Clut^ 1 ' S ) D0 ds +K9 f 0 C W S1 - 8 > D0 43 + 


K llJ^ c L at f s l- s ) J ^ 9 48 - K 12 °L 6a ( s l " S ) I)20 ds + 

p s i p si 

Ki 3 J Cl^Csi-sJDS ds + Ki4 Ql^C si - s)I> 25 ds (5) 


where the coefficients of the integrals are 


dC T 


„ „ \&x It St 

K2 Kl ( dC L \ 71 S 

yda /■«• 
d€i 

K 3 = " K2 da 1 

/(1Ql\ 

Ida /t St 
K 5 = Ki ^ /dQL\ 11 S 
\da / -y 

cUn 

*6 = “ K 5 £T 

'dCL \ 

, da /t St 
^ = K T /dCL\ ^ S 
i^da 


(5a) 


(5b) 


(5c) 


(5d) 


(5e) 
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Kg = -Kg 


del 

da 


(5f) 


Kll 


dCL\ 'St P 2 U 
da /t- 11 W 2 V c 


(5g) 


Kis 


/ \ St P y 2 Zy del 

Ida / ^ W 2 c da 


(5h) 


k 13 


>( 3.Cl \ p St v 2 toe 
to jt 2 W dS 


(5i) 


Ki4 


£ St v 2 £t At- /SCl\ t 
2 ¥ c At~+ 2 \§D8 / a 


'dCL\ /SCl\ 

^ jt \^8 /b 


(5J) 


to e 

Either experimental or theoretical values may he used for 

te,\ d6 

appearing in equation (5i) and r^gg). and [gggj In equation ( 5 J ) - 

’ A \ B 


Some theoretical values are given in figure 1(a) of reference 3. Other 
expressions appearing in equations (5) and (pa) to (5,3) are derived in 
the section entitled "Aerodynamic Coefficients." With the load-factor 
increments An w and Ant expressed as a function of-the aerodynamic 


coefficients and airplane configuration, equation (2) can he expressed 
in terms of the airplane parameters. All hody moments are taken about 
the center of gravity and assumed positive in the nose-up direction. 


The hinge moment due to the aerodynamic forces, assumed positive 
when the trailing edge is depressed, is set equal to the inertia hinge 
moment- taken about the axis of rotation of the elevator. An expression 
for the total hinge moment H in terms of- the aerodynamic character- 
istics and configuration of the airplane is required. The total hinge 
moment is made up of the aerodynamic hinge moment due to a change in 
angle of attack on the tail and the hinge moment - due to elevator 
deflection. 
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The equation for the total hinge moment is as follows: 


U 


' S 1 


H = " Kl 5 + K l6 J An ds - K 1T a + K l8 D9 + 


K19 


r B 1 p B i 

/ e a (si - s)An(s) ds - Kgo / e a ( B l - s)D9 ds + 
J 0 J 0 


r s i r s i 

K21 / €q(si - s)D^9 ds - K22 / sg(si - 


U 


s)D - ds + 

Umax 


K236 - K24 D 5 

where the coefficients are 


( 6 ) 


K15 - - 2 UmaxVSeQt^ce (6a) 

Kl6 = - | cSeCegCho, (6b) 

KlT - - | V 2 S e c e Cha (6c) 

K18 = | V 2 S e c e Ch ct ^ (6d; 

p 

K19 = 2 cSeCeCi^ ^ g (6e) 

K20 = | V 2 S e c e Ch a ^ (6f) 

%1 - - I V^eCetV || h (fig) 

^22 = ^ UmaxVSeceCi^ ^ (fill) 
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C 


K24 = % v2 Se c e 


At 


A t + 2 


K2 3 = | V 2 S e c e C hs 


/^h\ 

\^D5 /a \da< /t\^SS /b 


ft + 4 q 2 k 

c pVS e CeC 


(61) 

<6J) 


The Duhamel's integrals appearing in equation. (6) are not due to 
any unsteady hinge moments, hut merely afford a means of- expressing 
the effects due to downwash. These integrals establish a relationship 
between the changes of~angle of attack on the tail caused by the changes 
of angle of attack on the wing. If the value of the expression in 
brackets in equation (6j) cannot be determined experimentally, its 
value is given theoretically with the aid of figure 1(b) of reference 3» 


For the purpose of simplifying the differential equations. 


equations (2) and (3) are rewritten as follows: 

D 2 0 = K25 An w - K£6 Ant "'j 

K27 iPd +-K28 3^8 + K29 An = H ! 

> (T) 

where the coefficients are 

m^c 2 
Ke5 ■ IV 2 

(7a) 

Wl + c 2 
K26= IV 2 

(7b) 

K27 = j 

( \ y2 

% + HeZh - rI 8 U- 
^ ' c ^ 

(7c) 

k 28 = | 

'ie + r 2 !^ — 
S / C 2 

(7d) 

Kg 9 = • 

-(Hg + rH e )g 

(7®) 
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Aerodynamic Coefficients 

Unsteady lift .- Theoretical unsteady-lift functions due to both a 
unit jump of angle of attack and a penetration of a sharp-edge gust are 
given in reference 1 for elliptical surfaces of several aspect ratios 
for incompressible flow. The normalized functions are reproduced here 
in terms of the whole reference chord of the wing as distinguished from 
the half -chord notation used in reference 1. 

For a unit jump of angle of attack, 


Q L(X (s) a _ 3 = 1.000 0,283e rl * 08 ° s 

Ct(b) a c = 1.000 - 0.36le"°' T62s 

CL A=o 



C La (s) A=00 = 1-000 - O.^e" 0,0903 - 0.335e‘°' 6 ° OS 


J 


For a penetration of a sharp -edge gust, 
pL g (s) A=3 - 1,000 - 0,679e" 1 - ll6s - 0.227e” 6,liOs 

QL g (s) A= g = 1.000 - 0.4k8e"°’580 s - o.272e“ 1 ' 1< -50 s 
C L s ( s )a=«- = 1 * 000 ' 0.236e"°' 1:L6s - 0.513e^°' 728s 




- 0.193e" 6 '°° s 

-4.84s 


> 


0.171* 


J 


(9) 


The expressions for the finite -aspect-ratio unsteady-lift functions 
were derived using the midchord of the elliptical surface as the 
reference chord. 


In establishing the unsteady-lift functions on the tail, the 
exponents of the foregoing expressions must be multiplied by the 
ratio c/c-fc in order to convert these expressions from tail chords to 
wing chords since they involve the distance traveled by the airplane 
measured In terms of the wing chord. 


Unsteady downwash . - The theoretical value of the Uh st eady - dovnva sh 
angle at the horizontal tail due tg a unit jump of angle of attack on 
the wing is obtained by evaluating equation (1) of reference 4 for an 



Ik 


KACA TN 2035 


average configuration. The result may be normalized and approximated 
by step functions to give the desired function e a . However, the curve 
in figure 2(a) is seen to be principally a function of the distance 
between the trailing edge of the wing and the leading edge of" the tail 
and this distance, due to configuration plus an additional aerodynamic 
lag, account's - for the interval between negative and positive step' 
functions. For the unsteady -downwash-angle function due to penetration 
of a gust. €g, an. additional lag of 0.75 chord is applied to these 
step functions. (See fig. 2(b).) Figure 9 of reference 4 indicates 
that tail length does not affect the magnitude of the function appreci- 
ably butr merely tends to affect the interval of discontinuity. In 
order to generalize the functions e a and eg with respect to air- 
plane-tail length, the lag interval is considered as a sum of the 
distance from the trailing edge of the wing to the leading edge of the 
^ tail plus an additional aerodynamic lag of 0.23 chord. The steady- 
state values of ea and eg were found to be approximately equal to 
the value for de/da. 

The effective unsteady-lift function on the tail due to the down- 
wash from a unit jump of angle of attack on the wing may be determined 
by a similar procedure. With the results obtained from the solution 
of equation (1) of reference ij-,- equation (2) in this same reference can 
be evaluated. The results' for an average configuration, when normalized 
and approximated by step functions in a manner referred to previously, 
give the desired function QLe a - (See fig. 3(a).) The effective 

unsteady-lift function due to downwash which results from the penetra- 
tion of the wing in a sharp-edge gust Cl € is obtained by applying an 

additional lag of 0.75 chord to the step functions describing C;L 6a . 

(See fig. 3(b).) By an analysis similar to the one used to generalize 
the functions for unsteady -downwash angles, these two unsteady-lift 
functions are generalized with respect to tail length by considering 
the total -lag interval from negative to positive step functions as a 
sum of the distance fronuthe trailing edge of the wing to the leading 
edge. of - the tail, plus an additional aerodynamic lag of 0.80 chord. 

The steady-state value deq/da is approximately equal to 90 percent of 
the value for de/da. 


Transformation to Matrices 

Stick -free stability .- The integral relationships required for 
transforming the differential equation into a form for numerical solution 
are given in appendix A. In the process of obtaining the numerical 
solution, which is in the form of a recurrence formula, a number of 
constants were used to simplify the writing of the final solution. 

These constants, which appear as elements of a square matrix and several 
rectangular matrices, are in terms of the airplane stability 
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coefficients Ki to K29, the transformed unsteady-lift functions B 
and C, and the increment As. The constants Kq to K29 are defined 
in appendix B. The constants given in appendix B are for a six-point 
solution, that is, for a recurrence formula containing six previously 
known values of the variables to be solved. A six-point solution was 
selected in this paper since itj is believed that the transformed 
unsteady -lift functions B and C can, be approximated by six values 
at 1.5-chord intervals with sufficient accuracy even though only a 
moderate expenditure of time is required. The accuracy of the numerical 
solution would largely depend upon the size of the interval and the 
rapidity of the motions that the airplane undergoes when disturbed from 
equilibrium. For airplanes with a high mass parameter, a high moment 
of inertia about the center of gravity, and a high moment of inertia 
for the elevator and stick, the accuracy of the method is expected to 
be improved for the given six-point solution. 


Certain terms appearing in the equations for the E 
stants defined in appendix B have the subscript m - Zcl 


and F con- 
and certain 


terms in the equations for the G constants in this same appendix have 


£wt 

c 


+ 0.80 


Zwt 

c 


+ 0.23 


As 


the subscript m - Z e , where Zql = — and Z 6 = 

The physical significance of these terms in the equations for the 
constants is attributed to a lag of tail forces and hinge moments 
existing on the tail which were caused by the effects of the unsteady- 
downwash functions. (See figs. 2(a) and 3(a).) In order to express 
these effects mathematically, two values of the constant must be 
evaluated, one value with the term containing the subscript neglected 
and the other value with this term retained in its calculation. For 
values of m < and m <Z € , the first values of the constants are 
used in the calculations} for extended values of m, that is, m ^ Zcl 
and m = Z e , the second value is used. Since the numerical evaluation 
is made for a predetermined interval, the amount of lag specified by 
the discrepancies in the unsteady-downvash functions cannot conveniently 
be taken into account exactly} therefore, the amount of lag is taken 
to the nearest value of m. 


When the constants E, F, and G are assembled into a square 
matrix and a series of rectangular matrices, the three transformed 
equations may be solved s imult aneously by an inversion of the square 
matrix. The variables describing the motion of the airplane can then 
be solved for the stick-free condition by the following recurrence 
formula: 
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r 


fAn ^ 


(d 2 ©/ 

p2s 


m 


U 1 

u 3j 


V 


+ s 


m 


An ' 

d 2 © 

jpb 

•jAn 

jD 2 © 

t d2& 

r D 2 © 


ii n 
ii® 28 


y m-1 


An 

D 2 © 

Jpb l 

-j-d 2 © 

r D^8 


' m-2 


\ 


(An ^ 
D 2 © 
D^s , 
j-D 2 © 

r D ^5 


> do) 


ni' 


-6 




where the forcing-function matrix ( 


Ull 

Uq? is defined in appendix C for 


U' 


3J 


either a sharp-edge gust or. one of arbitrary shape. With the aid of 
equations derived in appendix A, equation (10) is evaluated in a manner 
analogous to the solution of the wing penetrating a sharp-edge gust 
described in the section of the paper entitled "Solution of One Degree 
of Freedom." 


The rectangular matrices QeQ are evaluated together with the 

square matrix [m] in terms of the constants E, F, and G defined in 
appendix B. The square matrix represents the inverse of the matrix of 
the coefficients of An, D 2 ©, ana D 2 6 at the point m. 

Stick-fixed stability . - The motions of the airplane with the 
elevator- in a locked position may be determined from the stick-free 
stability equations by modifying the matrices appearing in equation (10). 
The effect of eliminating the elevator motions is to eliminate in the 
numerical solution certain E and F constants and all the G constants. 
Also, appropriate changes in the recurrence formula are made by 
modifying the square and rectangular matrices and the column matrices 
containing the variables and their integrations describing the elevator 
motions. In appendix B the modified values of the matrices Tm]o and CsIq 
are presented and are written with a subscript zero to differentiate 
them from the stick-free condition. 

The forcing functions Ui and U2 remain unchanged whereas the 
function U3 is eliminated. In accordance with the previously cited 
changes, the following recurrence formula is obtained and can be 
evaluated in a manner similar to that for the sick-free condition: 
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APPLICATION OF METHOD 
Solution of One Degree of Freedom. 


The comparison of a response to a dynamic system calculated exactly 
to the response calculated by the numerical method derived in appendix A 
will help clarify this method as well as give an estimate of the accuracy 
to be expected when applied to a more complicated system. In order to 
demonstrate the validity of the numerical integration method, a response 
to penetration of a sharp-edge gust has been calculated for a wing 
restrained in all but the vertical direction. The change in lift 
coefficient due to the gust and due to motion has been calculated as a 
function of wing penetration. Assuming a unit jump in angle of attack 

due to the gust ^ = 1 gives the equation of motion as 


dCr 1 ps l 

Cl(si) = - — j/^ QL a (si - s)Cl(s) ds (12) 


where p g is defined as the mass parameter of the wing as follows: 


= 


2W- 


aPL a 

g £T Spc 


In reference 1, Jones has solved equation (12) by operational 
methods for different -values of mass parameter p g . All calculations 
were made for a wing of aspect ratio 6. In reference 1 the unsteady- 
lift functions C L g and CL a are defined as functions of half-chord 

while in this paper they are a function of whole chords as well as 
being expressed as a fraction of the steady-state value. The exact 
solution of the response using the unsteady-lift functions for aspect 
ratio 6 and a value of p g = 13.2 (or in Jones* notation p = 124) is 
shown in figure 12 of reference 1. This curve is reproduced herein in 
figure 4 as a function of whole chords. 
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The transformation of -equation (12) to the following recurrence 
formula can be accomplished, with the aid of the equations found in 
appendix A: 


'PL* " 
CLm-l 
^^01-2 



° L m-5 

CL m-6 i 


where 

^ = Vl + As (fc + H C Lm-l “ 2k ^-2 + 

The following numerical values are substituted in equation (13): 


ac L 

da 


1.50n 


-i- = 0.0760 

(ig 

As = 0.75 chord 
B(s) - 1 - ClJb)^ 

(where the values Bg, Bp, - . . . Bg are computed for the inter- 
val As =0.75 chord from normalized values of the function QLo.( s )a=6 

shown in fig. 9 of reference 1). The reason for choosing a smaller 
interval is that the mass parameter pg is very low. The resulting 
recurrence formula is given as follows: 


- °-°75T6l - o.oa^CLfc-a + 0 - 01&( tn-2 + 

0.00256CL m _ 3 + 0.00l40CL m _ 1( + 0.0016701^ + 0.0C<)265CL m _gj (15) 
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The forcing function CLg^ is computed for an interval 

of As = 0.75 chord from normalized values of the curve of figure 9 of 
reference 1. Equation (15) can he ccmputed step by step in a manner 
similar to that described previously. The results of these numerical 
computations are shown in figure 4 together with the exact solution 
obtained by operational methods. 


Solution of Stick-Fixed Condition 

In order to establish a general procedure for performing the 
necessary calculations in a stability analysis, an illustrative example 
is presented in detail. The aerodynamic characteristics corresponding 
to the configuration of a modem transport airplane were selected. The 
parameters necessary for a stick-fixed analysis are shown in table I. 
Values of It and 2y axe given for three center -of -gravity positions. 
Accordingly, the stability constants Kq to K 29 may be evaluated by 
equations (4a) to (4d), ( 5 a) to (5j)., ( 6 a) to ( 6 j), and ( 7 a) to ( 7 e). 

The results are shown in table II for those constants involving only 
the stick-fixed condition. 

The unsteady-lift functions due to a unit jump of angle of attack 
on the wing or tail were calculated from equation ( 8 ). Note that the 
tail unsteady-lift functions must be expressed in wing chords by 
multiplying the value of s by c/c-t- The B(s) and C(s) functions 
for infinite aspect ratio for an interval of As = 1.5 chords is shown 
in table III, Together with the constants Kq to K 29 , B(s), C(s), 
and As, the E, F, and G elements can be calculated with the use of the 
equations in appendix B. These elements of the matrices QmJ and [sj 
are shown both for the stick-free condition and for the stick-fixed 
condition in this appendix. Accordingly, the values for the constants E 
and F required for the stick-fixed condition are shown in table IV for 
an airplane with the center of gravity at 25 percent of the wing mean 
aerodynamic chord. The E and F elements are given as a function of the 
variable m. in the recurrence formula. None of the G elements are 
required for the stick-fixed condition. Since a discontinuity due to the 
effective unsteady -lift coefficients on the tail caused by a unit jump 
of angle of attack on the wing CL eQ . occurs at m = 2 , two values of 

the elements Eq, Fq, and Fg appear in the matrix CmHq - Consequently, 
two values of the matrix DCq have to be calculated. 
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The forcing functions represented by the column matrix 



can be 


calculated by equation (Cl) for a sharp-edge gust or by equations (C2) 

and (C3) for an arbitrary- shape gust. The unsteady-lift functions due 

to penetration of a sharp-edge gust Cl~(s) were calculated by equa- 

' & 


tions (9). By a method similar to the one used for the unsteady-lift 
functions due to a unit jump of angle of-attack on the wing, the tail 
unsteady-lift function due to penetration of a sharp-edge gust CL gt (s) 


is expressed in wing chords by multiplying the value of s by c/cq. 
The effective normalized unsteady-lift coefficient on the tail due to 
wing penetration Cl Pi ,(s) is shown in figure 3(h). Although the dis- 

continuity in this function for the configuration considered occurs 
at m = 3 to the nearest integer, the value at m = 2 was chosen for 
simplification since the length of the discontinuity determined by 
reference 4 is not very accurate and the length of discontinuity due to 


lag of tail penetration r = • ^ ■ is equal to 2 to the nearest 

integer. Therefore, all discontinuities are made to occur at m = 2 
in the sample calculations. The forcing functions Uq, and Ug have 
been calculated for a sharp-edge gust - for-an airplane with the center 
of gravity at 25 percent of the mean aerodynamic chord and the results 
are shown in table V. 


The responses were calculated by equation (ll) together with the 
formulas given by equations (A5) and (A8) of-appendix A in a manner 
previously outlined for the calculation of the simplified example. In 
figure 5 the time histories of acceleration increment and angle of 
pitch about the center of gravity of the airplane penetrating a sharp- 
edge gust are shown from an evaluation of the recurrence formula for 
an average center -of -gravity position. The load-factor increment- at - 
the center of gravity can be broken up into wing and tail load-factor 
increments by evaluating equation (4) for determining An w (s) and 
then subtracting that function from An(s) to obtain the load -factor 
increment on the tall Ant(s). The results of“these calculations are 
also shown in figure 5- 

Although sample calculations are not given for the two other 

center -of -gravity positions, the results are shown in the form of 
maximum load-factor increments An^x for a sharp-edge gust and two 
gradient flat-topped gusts at three center-of-gravity positions in 
figure 6. The responses to the gradient gusts can either be calculated 
by equations (C2) and (C3) or by the alternate method previously described, 
which in this case would be the more expeditious since three gust-shapes 
are considered in this paper. 
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DISCUSSION 


An estimate of the accuracy to he expected in the solution of the 
equations of motion by the matrix method may be determined by an inspec- 
tion of figure 4 where both the exact and numerical solutions are given 
for the lift coefficient of a wing penetrating a sharp-edge gust. The 
numerical solution approximates the exact solution very well, except for 
the discrepancies occurring at the beginning of the curve. This condi- 
tion will exist in most cases because the first point of the motion is 
computed from only those terms at time m of the recurrence equation, 
the second point is computed from those terms at times m and m - 1, 
and the third point is computed from those terms at times m, m - 1, 
and m - 2. In like manner, the other points can be computed until the 
full recurrence equation is in use. The error reduces as the mass 
parameter increases since the motion for the first few calculated values 
is reduced. Consequently, for low mass parameters where the motion is 
rapid a smaller increment of As than 1.5 chords must be taken. If 
smaller increments are required, the values of B(s) and C(s) may 
have to be approximated by more than just six points, otherwise the 
assumption that these functions are zero after the last point is not 
justified. In order to overcome this difficulty, the recurrence 
equation should be calculated from m = 0 to m = 6 for an increment 
of As less than 1.5 chords. The recurrence formula is then reevaluated 
for an increment of As = 1.5 chords by using values of the variables 
and their corresponding integrals taken from the first calculation. In 
addition to the effects of mass parameter, the accuracy is also a 
function, of the forcing function and their corresponding integrals 
taken from the first calculation at intervals of As = 1.5 chords. For 
forcing functions having a large initial slope, the motion for the first 
few steps would again be large and therefore the errors corresponding 
to its calculation would be substantially increased. As before, the 
errors may be minimi zed by using a larger number of small increments 
of As. Attention to the details of calculation for the sources of 
error pointed out, however, should lead to accuracies of calculation 
by the matrix method for the stick-fixed and stick-free conditions that 
are at least as good as those shown for the single-degree-of -freedom 
case in figure k. 

The assumptions made in the derivation of the differential equations 
impose certain limitations other than that of confining the treatment 
to longitudinal -stability effects alone. The length of the transient 
response to be calculated Is limited since the effect of variations in 
forward speed was neglected. The inclusion of variations in forward 
speed, however, would unduly complicate the problem at the present time. 
Another limitation to consider is the effect the pilot has on the 
handling of the airplane. If the pilot's reactions are proportional 
either to the accelerations, velocities, rate of change of acceleration, 
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or to various combinations of them, the restraint can be included in 
the calculation of the airplane motions by assuming these effectB on 
the stick are due to an automatic pilot-; However, knowledge of the 
pilot's behavior when flying through. .turbulent air is insufficient— at — 
the present time to warrant any such representation. 

Examination of the sample total -load-factor response curve in 
figure 5 shows a discontinuity occurring approximately at the time when 
the tail first enters the gust. As seen, the discontinuity arises from 
the load calculated to be on the tail and is a result of the combination 
of the downwash effects on the tail with the effects of-gust penetra- 
tion. Although this' discontinuity in the response does not actually 
exist, it is allowed to remain in this form for the determination of" 
the response to arbitrary gust shapes when the response to the sharp - 
edge gust is used as an indicial response. In this way, the new 
response does not depend upon the manner in which the discontinuity is 
faired. 

The curves shown in figure 6 illustrate the results of the use of 
the method for determining the effect - of center -of -gravity position and 
gust-gradient distance on the gust -load response of a particular air- 
plane. The curves show that the effect of- airplane stability may be 
important in determining the gust load factor on an airplane. The gust- 
load-factor increment can be appreciably affected by the center-of- 
gravity position when the airplane is traversing large gradient gusts. 

A rearward center -of -gravity position, representing a decrease in 
stability, increases the maximum load-factor increment Anmay. On the 
other hand, with the center of gravity in a forward position , the 
stability is increased and the corresponding maximum load-factor 
increment-1.6 decreased. 

All the computations involved in the solution of the complicated 
integral equations of motion of an airplane penetrating a gust have been 
reduced by the method of this paper to simple, but lengthy, arithmetic 
operations. The use of some type of high-speed automatic computing 
machine (such as the Bell Telephone Laboratories X- 66 7^4 relay computer 
In use at the Langley Laboratory) would help overcome the difficulty 
of performing these computations. If the problem is set up in a machine, 
the length of time required for the computation of the stick-fixed 
stability equations for a single condition would usually require some- ■ 
what less than 1 hour to approximately 8 hours, depending on the type 
of machine. The initial time required for setting up the problem on 
the machine would appear less significant when the number of-solutions 
of the equations Is increased. 
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CONCLUDING REMARKS 


The numerical method derived in matrix notation for solving the 
differential equations of motion provides a means for determining the 
effect of the various stability parameters on the gust load factor in 
a reasonable time considering the complexity of the problem. Although 
the method predicts the total load-factor increment on the airplane, 
separate wing and tail loads may be computed with comparative ease. 


Langley Aeronautical Laboratory 

National Advisory C ommi ttee for Aeronautics 

Langley Air Force Base, Va., November 8, 19^9 



2b 


NACA TN 2035 


APPENDIX A 

SOLUTION OF EQUATIONS BY A METHOD OF NUMERICAL INTEGRATIONS 


The solution of the integral equations may he accomplished, by the 
calculus of finite differences. A recurrence formula is obtained In 
this manner for the variables in terms of their preceding values and 
the integrals of the preceding values, 


Two mathematical operations are required for a solution of the 
equations by the method considered in this paper. Firsty a numerical 
evaluation of the various Duhamel's integrals expressed as functions 
of the variables describing the motions of - the airplane is required. 
Second, a method is necessary to evaluate numerically a relationship 
between the variables and their respective first and second derivatives. 


Since the indicial response Cl^s) appearing in the Duhamel's 

integrals approaches unity asymptotically, a convenient substitution can 
be employed. so that- the required integrations for succeeding intervals 
need ncjt be continued indefinitely. For this simplification the 
substitutions required for the unsteady-lift functions are: 


B(s) = 1 - Cl^s) 
C(s) = 1 - 0 Lo&fc (e) 


(Al) 


The functions B(s) and C(s) approach zero quite rapidly; therefore, 
the integration of the product of this function with the derivative of 
the arbitrary function may be carried out in a few steps. If the prod- 
uct of- the functions to be integrated be y(s), the desired value of 
the integral at point m can be evaluated In recurrence form for six 
intervals with the aid of Simpson's rule by the following equation: 


Ym 

Ym-1 



ds 


_ _y _ A 1^242411 

- % as |_3 , 3 , 3 , 3 , 3 , 3 , H 


Ym-2 
( Ym-3 
7m-b 
Ym-5 


y m -6 


(A2) 


where 


As is the interval of integration of the function y(s). 
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The Duhamel’ s integrals containing the downwash functions may he 
simplified, to ordinary type of integrals,, because step functions 
represent the indicial responses. The integrals containing these 
functions may be replaced by two simple integrals, differing only in 
the upper limits of integration by an amount specified in the dis- 
continuity occurring in the indicial responses. As an example, the 
following Duhamel 1 s integral may be replaced by two integrals (see 

fig. (3)): 



The evaluation of the Duhamel 's integral now merely reduces to the 
second problem, namely that of relating the variable to its first 
derivative since the integration of the load factor is the first 
derivative of the displacement divided by the acceleration of gravity. 

When approximating the successive derivatives of a function 
numerically, the accuracy of the operations decreases for each operation 
since the curve approximating each derivative would have to be of a 
lower order. As an example, if a function is replaced numerically by 
parabolic segments, the numerical approximations to the first derivative 
would be a series of straight-line segments; the second derivative, 
obviously, would be meaningless since its numerical approximation would 
be a constant for any value of the original function. On the other 
hand, a numerical method consisting of successive integrations will 
tend to reduce the error for each series of operations because the 
curve approximating each successive integration would be of a higher 
order. Thus, if a function is replaced by parabolic segments, the 
numerical approximation to the first integration would be a series of 
cubic segments and the numerical approximation to the second integration 
would be a. series of quart ic segments. Consequently, for each con- 
secutive integration, the accuracy of the numerical approximations tends 
to increase. 

In the numerical solution, the highest derivatives of the variables 
can be conveniently treated in the differential equations. Any 
derivatives of lower order may then be obtained by successive 
numerical integrations. In. the previously derived equations of motion, 
no derivatives of higher order than the second occur; therefore, two 
successive integrations are required to establish the relationship 
between the variables and their respective derivatives. 
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may. be 


m 


When the function and its integral from 0 to- %_! are known, the 
value of the integral from the increment %-! to sm may be obtained 
from the known values y m _-j_, y m _2> y m _ 3 > . . . , and the unknown 
value y m . The integration over the increment from s m _i to Sm can 
be perfprmed in- a manner analogous to that used in Simpson's rule. In 
lieu of passing a parabola through these points, a cubic assumed to 
pass through four points may be used to accomplish the same purpose 
with increased accuracy. Equation (a 4) may be rewritten 


The first numerical integration of the variable y(s) 
performed in the following manner:' * 


Sm 


y(s) ds = 


sm -1 


y(s) ds + 


0 


0 


I 


s m 


y(s) ds 


Sjn-l 



+ As 


9 19 .5 1 

2P2F ££’2? 


A 


Ym 

ym-i 

^m-2 

y m-3 




(A5) 


For the double integration, an analysis similar to the one 
previously derived may be utilized 



y(s) ds ds = 



I y ds 



I y ds + 


P * 


ds 


^-l 


(A 6) 


or, in terns of the double integration symbol II y , 



* Ag [^’H’-2T'5E 



m-1 

m-2 



(A7) 


If the value of from equation (A5) is substituted in equa- 

tion (A7), the double integration of the function y(s) can be 
written in terms of the previously determined values of the function, 
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its first and second integrations, and the function y(s) at the 
next desired value s^ as follows: 






r il + 


<As)‘ 


9 

9 19 5 l 


|ym-ll 


7 5 1 

2? 

jp2?'"25>25 

< 

|y m -2| 

,> + As 

5* “25*25 


y m -3 


'4-i' 

i-3 


(AS) 


The solution of the equations of motion by equations (A5) and (a 8), 
together with the evaluation of Duhamel's integral with the aid of 
equation (A2), can now be accomplished by a systematic procedure. 
Complete histories of the accelerations, velocities, and displacements 
of the airplane are computed simply by evaluating, respectively, the 
variables y m , 3^ and H^. The procedure can be continued for as 
many steps as desired. 


In order to satisfy the initial conditions of the differential 
equations, the initial values of displacements, velocities, and 
accelerations may be substituted into the numerical solutions. For the 
problem considered in this paper, all the initial values of the dis- 
placements and velocities are zero, since the airplane is assumed to 
be in steady level flight before being disturbed by the gust. If the 
gust velocity also increases from a zero value at time zero the initial 
accelerations are zero as well. The numerical solution may be regarded 
as a recurrence formula for the values y m , 3^, and 13^ which con- 


stitute the acceleration, velocity, and displacement, respectively. 

With the values at m together with those found at m - 1, m - 2, 
m - 3# • ■ • , the values at m + 1 can be found. In like manner, the 
other values of m can be obtained. 
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APPENDIX B 

CONSTANTS USED IN CALCULATIONS 


When the differential equations of motion are transformed for 
numerical evaluation, they can be written more simply by combining 
groups of constants due to the airplane parameters into one constant. 
These substitutions are made in equations (4) to (7) and are tabulated 
as follows for convenience of calculation: 


Kl - (?£fe) 
\da Jv 2 


P ys tynax 
W 


Kq = K ; 


M = 


acl 

,6a )t 

11 T 


•'(den 

\da / w 


<2 da 


% 




XT w ^ It S t 

K 5 " Kj + /dCy 11 S 
W / w 


dei 

Kg = -K 5 ^ 


K - (^k) § E v 2 
l \da / v W 2 V 


(OCl) 

v V W / 1 St 
K 8 " k 7 /cLCl\ n S 
\da j w 

d?i 

K 9 - -K8 st 
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%0 


6Ql\ S £ v 2 lx 
^da J v W 2 c 


K- 


11 


d£L\ St P 2 H 
da yt ^ W 2 c 


k 12 

/dCL\ St p 0 

”\da 

7t 71 w 2 

k 13 


p St „2 (ia e 

\daj 

,!?’ B 

Kl4 

P St 

P Ct At 

2 W 



c A t + 2 


c da 


(§) A + 




k 15 ~ u max VS e C h a c e 

Kl6 = cSeCegCt^ 

K 17 = V 2 S e c e Ch a 


Ki8 = | ^SeCeC^ 


Ki 9 = | cSeCeCha f| g 

%> = | g 

Kei = -| V^SeceC^ g £ 


K22 = | UmaxVSeCeChc f| 

*23 = | V^eCeCiig 


k 24 = | ^Sec e j 


At 


A t + 2 



f ctC L \ /9Ch\ 

W /tV^/B 


£t 4Q 2 K 
c pVS e c e c 
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K25 = 


WZ t c 2 
I V 2 


' \ 

+ " rI s/ ~2 

wi s )$ 


% = -CH e + rH g )g 

When the equations of motion are transformed to numerical form by- 
making use of the integral equations derived in appendix A , terms 
common to the same time interval, that is, at times m, m - 1, 
m - 2, . . . , are grouped into one constant. These constants are 
obtained for each time interval, each equation of motion, and each 
variable and its integrals describing the motion of- the airplane. The 
constants formed by this grouping procedure are elements in the square 
and rectangular matrices in the recurrence formula. 

The constants are grouped with respect to the three equations of 
motion. The following constants are due to the transformation of 
equation (l) into recurrence form: 

Ei = 1 + ^ As K4 - 0.158K 6 + (l.l58K6) m _ 2cL + K5 - | As (KI^Bq + K 5C0) 
E2 = AB ) 2 _K T " K8 + O.I58K9 - (l.l58K 9 ) m _2 CL . + |- As K T As ^ + 


Kg As ^ + K 10 - K X1 - 0.1 5 8K 12 + (l.l58K 12 L 2c l + 
j As (-K 10 B 0 + KjjCqJ 

e 3 ■ k Ab (- k 13 k Ab * k 13 T ■ K u) + K lA AB t 
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E k = -K 4 - K 5 + 0.158K6 - (l.l58K6j m _ lCL _i 

E5 = l As K T + K8 - O.I58K9 + (l. 158X9} m _ IcL _^] - | As (k t Bq + 4 ^ + 
K 8 C 0 + ^^l) ■ K 10 + K 11 + °- 158X12 - (l‘158K 12 } m _2 CL _i 
E6 = ^ As^K 7 ^ K8 + O.I58K9 - (1.158X9)^2^-21 - |as(x T B 2 + K8C2) 
E ? = ^ AsjK 7 + K0 - O.I58X9 +(1.158X9)^2^.^ - | As(k 7 B 3 + K3C3) 

38 = | As (~K 7 B 4 - K8Cit) 

’9 = ^ As (-XyB^ - K3C5) 

10 = \ As(-k 7 b 6 - x 8 c 6 ) 

11 = X 7 + Kq - 0.158Kg + (1.158X9)^2^^ 

12 = K i3 ■ y ■ | c i) + k i4 

13 = K i3 **(-h " I c s) 

s l4 = k 13 As (^i: " f c 3) 
e 15 = ' k 13 Ab I c 4 
E l 6 = ~ k 13 As ^ C 5 

e 1T = ~ k 13 As \ c 6 _ 

e 18 = H As -Kk ~ *5 + 0-158X5 - (l.l58K5) m _2 CL _i + | As(k 4 Bi + X5CJ 
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E 19 

E20 
e 21 : 
E 22 : 

E 23 1 
E 24- : 

E 2p 3 

E 2 6 = 

E 2J = 


h As K 4 + % - °- 1 5 8k 6 + ( 1 - 1 58K^ m _ 2ct _ 2 + | As^Bg + KjCg) 

gif As ^ " K 5 + 0-158K6 - ( 1 * 1 58K6) m _ Zc ^_ 3 + |as(k4B 3 + K5C3) 

^ As^Ki(B4 

I as(kp 5 + K5C5) 

\ As(K 4 B 6 + K^Cg) 


19 

2¥ 


(As) 2 K ? Ij- 
( 1.I58K9 2lf) m _2 c _i 


K 7 I + ^ lr ^ y - 0.158K9 |j- + 


■l li 


+ As 


-K 


10 gy + K io 3 ^ + K n gy - K n f c i + 


O.lSSK^ || - (l.l58K le g) 


m-Z CL -l 


^) 2 


2 K T h + k T y - *6 h + K 8 T + °- 158k 9 24 


2 _ _ 


(l.!5 8K 9 |j) m _ 


^Ct - 2 


5 


+ As jKio gif + K 10 3 B2 - Kn ^ - Kn | C 2 - 


0. 158x^2 ^ ^ 1 . 158X22 

^-Ky^ + KB^-KBy - 0 . 158X9 |j- + 


ZZ/n-iCb-z 

l / * -\2[^l ^9_ ®o , T ,^ _9^ Co ,,ro T ,_ 9 


• Z Cl ‘ 3 


+ As 


( 1 ' 158k 9 &) m . 

0.2.58x12 g . ( 1 . 158K12 if) m .i CL _ 3 
| As ( K 10 B i|- - »LL c 4) 


-KlO + K10 | B3 + Kn ^ - Kn | C3 + 


b 
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E 29 = j As (ki 0 B6 - KnC6) 

19/9 c 0 \ 4 

E 30 = 2? As K13 As - K13 As — + X1I4. 1 - X]^ As 3 Cl 

E 31 = 3+ ^ ( _K 13 As + K 13 ^ T ' KlIt ) “ Kl ^ ^ | C 2 

E 32 = ££ ( K 13 ^ 2K ' K 13 As y + K l0 ‘ ^ As J C 3 

E33 = -K^ As Cj^ 

lj. 

E^ij. = -K^lj. As ^ C5 

E 35 = " K l4 As ^ c 6 

E 36 = K 13 

The folio-wing constants axe due to the transformation of equation (2) 
into recurrence form: 

F 1 = ^ As[k 4 K 25 - k 5 k 26 + 0-158K5K26 - (l.l^X^ 6 ) m _ Zc J - 

j As ^25^0 - K5K26C0) 

F 2 = 1 + As ) -K 7 K 25 + k 8 k 26 “ 0.156X9X26 + ^1. 158X9X25 + 
If As [k 7 X2 5 As ^ + X 1Q X2 5 - XqK^ As y + X 11 X s6 + 0.158X 12 X 26 " 
(l.l58X 12 K 2 6) m _ ZcL + ^ As ^-X 10 X2 5 B 0 - K-qX^CoJ 

F 3 = 2% As ^13^6 As - X 13 X26 As — + X^j^) - K^Xge As y 

f 4 = -K4K25 + * 5 ^ 26 - 0-158X5X26 + (l.l 58 K 6 K 2 6 )m-l CT -l 

. Li 
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F5 = ^As|k T K 25 - K8K 2 6 + O.I58K9K26 - (l.l58K9K26) m _2 ct .^] + 

| As (-K^Kg^Bo - kKsjH^^ + KqK^^Cq + j - K 10 K^ 5 - K^Kgg - 

0 . 158 K 12 K 2 6 + (l.aja^Kgg)^^.! 

F 6 = |f As |Tk 7 K 2 5 + ^26 " 0.158K9K26 + (l.l58K 9 K 2 6)m-lc L -2] + 

| As (-K 7 K 25 B 2 + KQK^gCg) 

F 7 = ^As|K 7 K25 " ^26 + O.I58K9K26 - ^ 1 . 158K9K26) m- Z Cl - 3 - 
3 As(k 7 K 25 B 3 - KgKg^) 

F8 = | As (-K 7 K2 5 B4 + K8K26C4) 

F 9 = | As (-K-fK^B^ + KgKggC^j 

F 10 = 3 A® (" k 7 k 25 b 6 + k 8 k 26 c 6) 

F 11 - k 7 k 25 - W26 + 0.1581^6 - (i.i 58 K 9 Kggj m _j ( ^_ 1 
f 12 = “ K 1 3^6 As £ + * 13*26 ^ y + K 1 3^6 As | C 1 - * 14*26 

F 13 = As (k 13 K2 6 ^ + K 13 K2 6 I c 2 ) 

Fl4 = As (-K13K26 + K13K26 jC 3 ) 

f 15 = k 13 k 26 Ab § c 4 

p l6 = Ki 3 K 26 As J C9 
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f 17 = *13*26 As J c 6 

F18 = || As pKliK25 + K5K26 - O.I58K6K26 +.(l.l58K6Ks6)m- !c L -l] + 

| As^K^Bl “ K5K26C1) 

F 19 = A® {*4*25 + 0.1^8K^2£ - ( 1 . 158K5K26) m - 1 c L -2 “ *5*26] + 
|as(ki^ 25B2 - *5*26 c 2) 

F 20 = y As [*5*25 + *5*26 “ 0 *158KgK25 + (l . 158X5X25] m _ zcl- 3] + 

| As^X^Bg - K5K26C3) 
f 21 = f As (ki^B^ - K^^gCij. j 
f 22 = 5 A* 3 ( k 4 k 25 b 5 _ *5*2 6 c 5) 
f 23 = | As (*5*25 B 6 " *5*26 c 6) 

F 2k = H (As) 2 [*^*25 H “ k 7 k 25 y ■ k 8 k 2 6 25 + *8*26 y + O.158K9K25 

(1.158X9X26 |f) m _ lcL _^] + [*10*25 H + *10*25 I B X - K11K26 || + 

*11*26 | c l - 0.158X^X25 g + (l. 158x^5 g) 

F 25 = (As ) 2 J^Ky*25 + *7*25 ~J + *8*26 2 % '*8*26 y ■ 0.158X9X25 H 

(1.158X9X26 §z) m _ lc _ 2 ~1 + As [*10*25 ^ + *10*25 | *2 + *11*26 Jj; + 
2 5 / 5 \ I 

*11*26 3 C 2 + O.I58X22K25 ^ _ (1*158X12*26 2?/m-lpL-2 
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*26 = |j : (a B ) 2 [k 7 K 25 Ij - K 7 K2 5 ^ - X^g |j + XgK^g ^ + 0.156X3X26 If - 
(1-158X3X26 ^(m- 2,^-3 + As [^10X25 gif + X 10 X2 5 y Ej - KjjKgg gjf + 
*11*26 f c 3 - 0.158K 12 K26 + (1. 158X12*26 

F27 = | As ^KioKg^Bl*. + K11K26C4) 

F 28 = f As ( K 10 K 25 B 5 + ^i 1 K 2 6 c 5) 

F29 = ^ As^KioK 25B6 + K11K26C6) 

X30 - || As(-X 13 K26 As 9 . + K 13 K26 As ^ - K^g) + K^g As | C 7 

F 31 = As (x 13 %g AS If - KyKgg As y + K^^Kjg 'j + Kj^g As | Og 

F 32 - gif (-*13*26 As |f + As ^ - K^Kggj + K^g 4. | C3 

f 33 = K l4 K 26 ^ I c 4 

i». 

f 34 = K l^2.6 As 3 c 5 
1 

f 35 = ^^26 As ^ eg 
F36 = -K13K26 

The following constants are due to the transformation of equation (3) 
into recurrence form: 

G 1 = k 29 - jlfAs Kl 6 " °‘ 1 ^19 + (1*167X19) m _j^J 
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G2 = Kqj + ^ As^ 2 K 1? - O.I67K20 + (^1.167X20) m- 2 e 
As |%g - O.I67K21 + (^-*^7K2i)ni._2gj 
G3 = K28 “ K23 (as ^J 2 - Kfilj. As ^ 

C4 = K l6 " G • 167^19 + ( 1 • l 67 K]_^j m- 2 g -1 

^ Asjj-K^Y + O.I 67 K 20 ~ (l* 1^7^20^111- 2 g -l + %8 

(l.l 67 K 2 l) m - 2 6 -l 

Gg = Jjj- As jl^.7 - O.I67K20 + (1.167X20) m- 2 6 -2 
G7 = ^ As P 1T + O.I67K20 - ( I.167K20) m- 2 e - 3 


- 6 . 167X21 


7 


Gq = As g + K24 


G9 = -K23 As ^ 


G^q = K2 3 As 

°11 = || As j %6 " 0-l67K 19 + (l.l67K 19 ) a . l6 .Ij 

G -| p = As j-K]_g + 0 . l6TK 19 - ( 1 . l67K 19 ) _gj 

g 13 = ^ As {^16 “ 0.l67K i9 + (l . l67K 19 ] m _ 6 .3] 

G lk = H As [fl7 As fe + K l8 + 0.167X20 As 2% - (1.I67K20 A s |f)m-2 6 -l 
0.167X21 + (l.l 67 X 2 i ) m . Ze _3 
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°15 = 2k AS ^-7 As 2% ' K l8 " 0 - l6 7%j As ^ + ( 1 * l67K 20 As ^) m -Z e -2 + 

0 . l67Kgi - ( 1 . 167K21) m- Z 6 -2} 


°l6 = As [~%7 As ^ + K l8 + 0 . 167 %) As 

0 . 167K21 + ^1.167X2^ m _j € .3 


19 


J 17 = 2? As ( K 23 As 2% + ^ 

G l 8 ■ A 6 ( -k 23 As H- ■ ^ 

G 19 = As ^23 As + %] 

G 20 = K 23 

C21- " k 17 + O.167K20 " (l.l 67 K 2 o)m-Z e -l 


(1.167%) ^ |r)m-Z € -2 


The elements due to the constants E, F ; and G are put Into 
matrices which modify the variables and the integrals describing the 
motion of the airplane. The matrix describing the variables at time m 
iB a square matrix that must be inverted in order that the variables 
at- time m can be solved simultaneously. The elements for this 
inverted matrix form the square matrix [M] which appears in the 
recurrence equation (10) . The elements of the matrix QQ i- s defined 
in terms of the constants E, I, and G as follows: 


M-l Mg 




Mg 


Mrj Mg M<J 
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where 


Ml = 1,253 

" 02=3 
H 

E-i Go - GoEo 
N 3,3 

M2- 5223 

- E2G3 
IT 

FlE3 - ElF3 

m 3 = 

- F2E3 
N 

F1G2 - G1F2 
**7 " u 

G3F3 

-F1G3 

H 

E 2 Gi - E-jGg 
“8 = H 


EiF 2 

- F]_E 2 


M<j — 

N 


where 


N — G 2 E 3 + G^EgF ^ — G^FgE^ ~ =l02=3 — F^EgG^ 

The matrices describing the variables and the integrals at 
times m - 1 to m - 6 form a set of rectangular matrices. These 
rectangular matrices consist of the elements defined by the 
constants E, F, and G. The matrices are given for each time interval 
as follows: 



=18 

=24 

=30 

Elj. 

E 5 

=12 

=11 

=36 

M- 

=18 

=24 

=30 

=4 

f 5 

=12 ' 

=11 

=36 


^11 

G l4 

°17 

04 

°5 

°8 

G 21 

020 


e 19 

E 25 

e 31 

e 6 

e 13 




M- 

F 19 

=25 

*31 

=6 

P 13 





ii2 

G 15 

°18 

g 6 

G 9 _ 
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to 




E20 

%6 

E32 

Ey 

E lk 


= 

F 

*20 

. ^ 

F 32 

f t 

F 14 



> 

G l6 

> 

°7 

G 10_ 



Egi 

E2Y 

e 33 

E 8 

e 15 

M 

= 

*21 

F 2T 

F 33 

f 8 

F 15 



0 

0 

0 

0 

0 



E 22 

^28 

E 3 ^ 


e i6" 

K 


f 22 

F 2 8 

F 3 k 

f 9 

F l6 



0 

0 

0 

0 

0 



®23 

e 29 

e 35 

: E io 

e 17 

K 

a 

f 23 

F 29 

f 35 

F 10 

f 17 



0 

0 

0 

0 

0 


All the square and rectangular matrices which have been defined 
previously are for the stick -free condition. The recurrence formula 
given by equation (ll), however, is a solution for the stick-fixed 
conditions and can be obtained by modifying the stick-free matrices. 
This modification consists of^eliminating all the G constants and 
those E and F constants' in the rectangular matrices which are post- 
multiplied by the variables in the column matrices that represent the 
acceleration, velocity, and displacement of the elevator. The matrices 
for the stick-fixed condition are written with a subscript zero and are 
defined as follows: 
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f 2 



-E 2 



EjF g — 

EgFi 


g “ EgF 1 

BOo - 







-Fl 



E i 



E 1 F 2 — 

EgFi 


EjPg — E 2 F 1 



E l 8 

E2k 


E 5 En 


0 " 

_£l 8 

*2b 


f 5 Fh 



e 19 

E 25 

e 6 ~ 


[S 2 

0 = 

F 

LI19 

F 

25 

f 6 _ 




E 20 

E 26 

e? 


[%] 

0 = 

Z 20 

F 26 

F 7_ 




**21 

^7 

e 8 ~ 


N 

lo " 

E 21 . 

f 27 

e 8 J 




s 22 

e 28 

*9*1 


[S5_ 

Io = 

_£22 

F 28 

f 9_ 




e 23 

1 

e 29 

E 10 


[*: 

0 = 

, f 23 

F 29 

r i°j 
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APPENDIX C 
FORCING FUNCTIONS 


The forcing functions Up, U 2 , and U 3 are dependent upon the 
shape of the gust that the airplane is penetrating. If the unsteady- 
lift functions are known for the sharp -edge -gust condition, the 
evaluation of the forcing function is very simple and may he calculated 
for the sharp-edge -gust condition from the following equations: 


'ui 




' ' 


^ N 

u 2 

= < 

KlK 2 5CL gw 

> + 

-K 3 Ks6CL £g 


-K 2 K 26 CL gt - 

„ u 3. 

m 

. 0 > 

m 

, “ K 22% , 

m 

" k 15 

V / 


(Cl) 


where, in the recurrence formula, T = — ^ - represents the time 

interval that the airplane has to travel before the tail enters the 
gustr; For m < t, the terms having the subscript m - t are neglected 
in the calculation of the forcing function. On the other hand, 
for. m ^ t these terms are included. (The value of t must be taken 
to the nearest - integer. ) 


Because of the discontinuities of the functions and QL € g 

(see figs. 2 and 3 )> two values of the forcing function exist at both the 
07 ^ 0 75 

points m = Zrv + and m = l e + In order to overcome this 

As e As 

difficulty when evaluating the recurrence formula, two values of the 
forcing function are calculated at this point. The recurrence formula is 

evaluated at times m = lop + and m = lg + by use of the 

first value of the discontinuity and then reevaluated by use of the second 

value. Therefore, the recurrence formula gives two-valued solutions of 

' 0 75 0 75 

the equations of motion at the points m = Zcp, + - ^ s and m = Z c + - . 

In order to evaluate the formula at' time m + 1 and. at succeeding 

intervals, the first 'solution obtained from the recurrence formula is 

neglected, thereby the second one is used as the solution for succeeding 

calculations. The numerical solution to the equations of motion therefore 

contains discontinuities due to the functions Cg and CL e g. For most 

cases, however, the magnitude of l e and are approximately equal; 

therefore, when these two distances are taken to the nearest integer 

of m, i e can be taken equal to Z^ and o^J-y one discontinuity 

appears in the solution. 
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For an arbitrary gust shape , a number of Duhamel 1 s integrals have 
to be evaluated for determining the forcing functions. These operations 
may be performed graphically or by various numerical methods shown in 
reference 5. The corresponding forcing functions for the arbitrary gust 
shape may be calculated from the following equations: 


where 




U 2 > = 


Kl*w 


KlK25$w 


r 


u- 


m 


+ < 


m 


3 e w 


-K 3 K 26 ^e. 

-K2 2 5 


w 


> + 


m 


r t ' 


-K 2 K 26 ®t 


Umax 


U 


m-T 


U 


y *) -L 1 c v si 


> Umax „ 
s) — — — ds 


1 0L gt ( S i - >)- 

®<% (s) -f o 1 C Leg (Bl - S)- 

r s i d 

l(s) = / 6 (sp - s). 

•Jq & 


ds 


U 

Umax 

ds 


U 

Umax 

ds 


ds 


ds 


U 


Umax 

ds 


ds 




(C2) 


(C3) 


If the response of the airplane for a number of arbitrary gust 
shapes is required, an alternate method would be far more expeditious. 

By first calculating the motions of the airplane for the sharp-edge- 
gust condition, the motions for each degree of freedom may be considered 
as indicial responses in Duhamel* s integral. When the gust shape is 
utilized as the arbitrary forcing function, the response for each degree 
of freedom may be calculated by Duhamel' s integral in a man ne r 
analogous to that used in calculating Up, U 2 , and U 3 for the 
arbitrary gust shape. Although the calculation of the response for the 
sharp-edge -gust condition would require additional work, the extra work 
is negligible when two or more responses to arbitrary gust shapes sure 
desired. 
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TABLE I 

AIRPLANE PARAMETERS 



per radian 



radian 


S, sq ft 

St, sq ft 

p, slugs/cu ft 

V, ft/sec - 

4 

dgp 

dot 

Aerodynamic center for wing and tail, percent M.A.C 




ft: 


For c.g. 
For c.g. 
For c.g. 


at 12^ percent M.A.C 
at 25 percent M.A.C. 
at 37§ percent M.A.C 




ft: • 


For c.g. 
For c.g. 
For c.g. 


at 122 percent M.A.C 
at 25 percent M.A.C. 
at 37^ percent M.A.C 


l vt> « * * 
U, ft/sec . 

g, ft/sec^ 

I, slug-ft l 2 


c, ft 
c t , ft 
W, It 


5.56 


. 3-21 

■ 738 

. 275 

0.002378 
308 
. 0.85 

0.500 

25 


-1.26 


0 

1.26 


. 36.64 
* 35.38 
. 34.12 

. 25.98 
50 

. 32.2 
209,600 
10 

. 7-58 

38,000 
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TABLE II 

CALCULATED VALUES OF AIRPLANE STABILITY CONSTANTS INVOLVING 
ONLY THE STICK-FIXED CONDITION FOR CENTER OF GRAVITY 
AT 25 PERCENT M.A.C. 


Kl = 1.9772 

K8 = 2.2272 

K£ = 0.36156 

- K9 = -1.1136 

k 3 = -0.18078 

0 

11 

0 

= 0.041342 

% = 7.8798 

K5 = 0.0075599 

e 

11 

0 

K£ = -0.0037799 

& 

VJ1 

11 

0 

Ky = 12.180 

Kgg = O.OO67616 
\naca^ 
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TABLE IH 

TRANSFORMED UNSTEADI-LIFT FUNCTIONS DUE TO UNIT JUMP OF 
ANGIE OF ATTACK ON WING AND TAIL FOR INFINITE 
ASPECT RATIO FOR INTERVAL As = 1.5 CHORDS 


B(s) = 1 - CL^Cs) 

C(s) 

11 

H 

/ 

CO 

Bq = 0.500 

Co 

= 0.500 

b-l = 0.280 

C 1 

= 0.240 

b 2 = 0.181 

°2 

= 0.147 

b 3 = 0.133 

c 3 

= 0.106 

= 0.105 

c 4 

= 0.084 

b 5 = 0.088 

c 5 

= 0.069 

b 6 = 0.075 

c 6 

= 0.057 



TABLE XV 


YALOBB OF I AHD F KLBMEBT8 HT BASBICBB JXb ABD t]S]o *3 A FWICTIO* OF VABIABI^ m 


7CB AS URPLABE VITH CEBU® OF GRAYXTI AT £5 FKKUKBT M.A.C. 


I and F 

E - 1 

m - 2 

» - a“ 

■ - 3 

m - 4 

1 

■ 

B1 

1.0156 

1.0156 

1.0131 

1.0131 

1.0131 

1.0131 

n 

-.000010246 

-.ooocoaa46 

-.0000015975 

■ -.0000015973 

-.0000015975 

-.0000015975 

Ee 

-5-0506 

-5.0506 

-4.6426 

-4.6426 

-4.6426 

-4.64a6 


I.OI9674 

1.019674 

la 016916 

-.019499 

I.OI6916 

1.016916 

1.016916 


-.049499 

-.049499 

-.045122 

-.043122 

-.043122 

Ft 

.000055155 

.000055155 

.000055155 

.000025559 

.000025339 

.000025539 

*5 

£1.908 

21.908 

21,908 

19.651 

19.651 

I9.ea 

*3 

-.070717 

-.070717 

-.070717 

-059458 

-.035458 

-.035458 

*6 

-7.0683 

-7.0883 

-7.0883 

-7.0083 

-6.6853 

-6-6833 

r 6 

.0072836 

.OO7E056 

.0072856 

.OOTS856 

.0045608 

.0045608 

*7 

-£.8010 

-2.8310 

-2.8010 

-2.8010 

-2.8010 

■ -2.8816 

? 7 

.0081800 

.0021800 

.0021800 

.0021000 

.0021800 

.002725 


-1.4633 

-1.4655 

-i.4655 

-1.4635 

-1.4635 

-1.4655 

*e 

.O0126SO 

.0012620 

.0012620 

.0012620 

.0012620 

.0012620 . 

i 

- 2 . kA-yi 

-£.4493 

-2.4493 

-2.4493 

-2.4493 

-2.4493 

.0000662 

.0000662 

.0020662 

.0020662 

.0020662 

.0020662 

--5G083 

-.52025 

-.52025 

-.32023 

-.52005 

-.32025 


.000858390 

.000429193 

.000429195 

.000429193 

.000429195 

.ooo4e90-95 

14.583 

14-383 

14.583 

13.294 

13.294 

• 13.094 

FU 

-.Ol6e49 

-.016249 

-.016249 

-.0075297 

-.0073297 

-0075297 

*18 

-.0319969 

-.0319969 

-.0319969 

-.02 6799 

-.026799 

-.026799 

*18 

.000040940 

.000040940 

.00004094a 

.0000057943 

.0000037945 

.0000057943 


■ v£!»<?60 



mLa/a 

fiPp^Qp 


*1? 

-.000024730 

-.000024730 

-.000024730 

-.000024730 

-.000015481 

-.00002.5481 

*so 

.0095075 

.0093075 

.0093075 

.0095073 

.0093075 

.0097801 


-.0000073996 

.noLorkk 

-.0000073998 

^nnkovkl 

-.0000073998 

_ nrtkyrkS 

-.0000073998 

-.0000073996 

ooV^VV 

-.0006092496 

^Anlifi'74 L 

ffcl 

-.000004SB36 

-.0006042836 

-.00000428361 ■ 

-.0000042836 

-.0000042836; 

V -.0000042036 11 

*22 

.0083134 

.0083134 

.0083134 1 

.0083134 

.0083134 

.0083134 

?£g 

-.0000070133 

-.0000070133 

-.0000070133 

-.0000070133 

-.0000070133 

-.0000070133 


.0017658 

.0017658 

.0017658 

nnrrtWt 

.00176*58 


S 

-.00060145683 

-.000001436® 

-.00000149683 

-OOOOOI45685 

-.00000145603 

-.00000145683 

11.036 

11.036 

11.036 

10.174 

10.174 

10.174 

*24 

-.044077 

-.044037 

-.044057 

-.038233 

-.038233 

-.038233 ■ 

*£3 

-5-05S 


-3-0535 

-5-0355 

-4.8e80 

-4.8288 

*2*5 

.0261^-1 

•0C6l4l 

.026141 

.024608 

.024606 


-.89033 

.0076401 

-.89203 

-.89203 

-.89203 

-.89203 

-93737 

*26 

.0076401 

.0076401 

.0076401 

.0076401 

.0073333 

%7 

-.66033 

-.66033 

-66033 

-.66033 

-.66033 

-.66033 

F £7 

.0044649 

.0044649 

.0044649 

.0044649 

.0044649 

.0044649 

$ 

-.10811 

-.10811. 

-10811 

-.10811! 

-.10811 

-108U 


.0073100 

,0073100 

.0073100 

.0073160 

.0073100 

.0073100 

*09 - 

-.224575 

-.224573 

-.224375 

-.224575 

-.224375 

-.224575 

?» .1 

.OQ15105 

.0015135 

.0015105 

.0015W 

.0015165 

.OQI5IH9 


“Two ™iu*« occur «t b “ 8 bccftuu of the effect of the function Cl^ cb the tell, 

°Fbr ill amine a of * greeter than 5, nee the veluaa of E and T te the coluen for m - 3. 





i 


i 


NACA TS 2035 



TABLE V 


I « 


CALCULATED FORCING FUNCTIONS Ul AND U2 FOR A SHARP-EDGE GUST 
FOR AN ADR PLANE WITH CENTER OF GRAVITY AT 25 FERCEOT M.A.C. 


© 

© 

© 

® 

© 

® 

© 

<D 

(D 

m 

K1CL„ 

K2 CLg* 

ti-t 


(Ul)m -©+©+(£) 

K25© 

-K26® 

-K26® 

(U £ ) m -©+© + (§) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1.2456 

0 

.02856 

1.2742 

0 

0 

-.00019311 

-.00019311 

2 

1.5323 

0 

.02856 

1.5609 

. 0 

0 

-.00019311 

-.00019311 

2 

1.5323 

0 

-.18078 

1.3515 

0 

0 

.0012224 

.0012224 

3 

1.6628 

.24980 

-,18078 

1.7318 

0 

-.0016890 

.0012224 

-.0004666 


1.7301 

.29724 

-.18078 

1.8466 

0 

-.0020098 

.0012224 

-.0007874 

5 

1.7775 

.31626 

-.18078 

1.9130 

0 

-.0021384 

.0012224 

-.0009160 

6 

1.8111 

.32692 

-.18078 

1.9572 

0 

-.0022105 

.0012224 

-.0009881 

7 

1.8388 

.33433 

-.18078 

1.9924 

0 

-.0022606 

.0012224 

-.0010382 

8 

1.8605 

■34001 

-.18078 

2.0197 

0 

-.0022990 

.0012224 

-.0010766 

9 

I.8803 

.34442 

-.18078 

2.0439 

0 

-.0023288 

.0012224 

-.0011064 

10 

1.8922 

•34797 

-.18078 

2.059^ 

0 

-.0023528 

.0012224 

-.0011304 

li 

1.9040 

.35075 

-.18078 

2.0740 

0 

-.0023716 

.0012224 

-.0011492 

12 

1.9129 

.35295 

-.18078 

2.0851 

0 

-.0023865 

.0012224 

-.0011641 

13 

1.9219 

.35473 

-.18078 

2.0959 

0 

-.0023985 

.0012224 

-.0011761 

14 

1.9278 

. 35614 

-.18078 

2.1032 

0 

-.0024081 

.0012224 

-.0011857 

15 

1.9327 

.35726 

-.18078 

2.1092 

0 

-.0024156 

.0012224 

-.0011932 
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Figure /.- Forces and moments assumed oct/ng on 

o/rpbne dur/nq penetration of a gust. (Distances 
forces, and moments shown m postive direction) 
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Q / 2 <3 4 5 e 


s, chords 

Co) Due fo un/fjump of ongJe of of fa ck or? wing. 



O / 2 o’ 4 d 6 


s, chords 

(b) Due fo penefrof/on of iv/ng ?nfo a sharp -edge gust 

figure 2 . - /Vor/vaf/zed uns. feaalg- down wash angle at 
hor/zon faf fa//. 
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s, chords 

Co) Due Co unit jump of angle of attack on w/ng. 



s / chords 

(h) Due fo penefraf/on of w/ng mio a sharp -edge gust. 


F/gure 3 . - Effective /Tar/nat/zed unsteady-/ iff function 
on the fa//. 
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F/gure #■- Comparison of numer/caf and exact -sduf/onc 
of the J/ff co&ff/ctenf an a w/ng penefraf/ng 
a sharp-edge gust ja^ — /3.Z . 



6, rctof/ans 


VJl 



.05 t— ’ 

o 

“ 06 
-JO 

%/s _ 

6 ? /O 15 £0 £5, 30 35 W 

(b) Angie-of- pitch displacement. S/ c ‘ lor ^ s 

F /gone 5 - Typ/ca/ responses for t/ie penetration of a sharp-edge gust with on 
average cenfer-of~ growfy /oosit/on. 


£ i » 
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O forward c.g 
□ /tear ward a g. 
A A verage c. g 


Gradient distance , chords 


O Sharp-edge gust 
□ /O-chord gradient 
A ftychord gradient 


a g,^ percent M A, C. 


figure 6 . - fffect of gust shape and cenfer-of -gravity 
position on gust toad factor . 
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